
read "encoder.mpl":
read "decoder.mpl":
read "channel.mpl":

# modify here for different parity-check matrices.
read "gallager.816.1A3.455.mpl":
#the matrix needs to be named H.

# generate parity-check matrix with dimension M * N
(M, N) := LinearAlgebra:-Dimension(H);
K := N - M;
t := columnWeight(H);

# generate matrix with dimension K * N
G := computeG(H);

# rate >= R, if H is linearly independent, then rate = R.
R := K/N;

# maximum iteration number
maxiter := 30;

# seed the randomizer with system clock.
Seed := randomize();

# generate a random sequence of bits of length K as the transmitted
# codewords.
genRandomMessage := proc()
global G, K;
    encodeLDPC([seq((rand() mod 2), i=1..K)], G);
end proc:

# verify the correctness of the generator matrix G computed above
verifyG := proc()
global H;
local msg;
    msg := genRandomMessage();
    evalb(add(i, i=LinearAlgebra[Modular]:-Multiply(2, H, Vector(msg))) = 0)
end proc;

# compute the number of error-bits.
countErrorBits := proc(trans, corrected)
global N;
local i, result;
    result := 0;
    for i from 1 to N do
        if (trans[i] <> corrected[i]) then
            result := result + 1;
        end if;
    od;
    result
end proc:

# simulate LDPC codes on binary symmetric channel.
simulateBSC := proc(p, blkno)
global H, M, N, K, t, R, G, maxiter;
local i, trans, recv, corrected, decodedResult, successflag;
local errbits, totalErrBits, totalErrorBlocks, totalIters;
local correctblkno, undetblkno, detblkno, stat;
    # initializing counters
    correctblkno := 0;
    undetblkno := 0;
    detblkno := 0;
    totalErrBits := 0;
    totalErrorBlocks := 0;
    totalIters := 0;

    # start decoding
    for i from 1 to blkno do
        print("Now sending block:");
        print(i);
        trans := genRandomMessage();
        recv   := bsc(trans, p);
        # a pair is returned (flag, decodedMessage)
        # currently using sum-product algorithm (not in log-domain)
        decodedResult :=  decodeBSC(recv, H, maxiter, p, 2);
        successflag := decodedResult[1];
        totalIters := totalIters + decodedResult[2];
        corrected := decodedResult[3];
        errbits := countErrorBits(trans, corrected);
        if errbits = 0 then
            if successflag then
                correctblkno := correctblkno + 1;
            else # this branch is not quite possible to happen unless a bug.
                error "warning, a bug happens";
            end if;
        else #if errbits > 0 then
            totalErrorBlocks := totalErrorBlocks + 1;
            totalErrBits := totalErrBits + errbits;
            if successflag then
                # then this is an undetected error
                undetblkno := undetblkno + 1;
            else
                detblkno := detblkno + 1;
            end if;
        end if;
    od;

    # compute statistics
    # bit error rates, block error rates, average iteration number
    print("Simulation statistics:");
    stat := [evalf(totalErrBits/(N*blkno)), evalf(totalErrorBlocks/blkno),  evalf(totalIters/blkno), correctblkno, undetblkno, detblkno];
    print(stat);
    stat
end proc:

#time(simulateBSC(0.074, 1000));
